1 Problem A - Euler-Maruyama Method

First we implement a discretization for the Black-Scholes-Merton (BSM) price dynamics with the given parameter set, with \(n = 250\) steps. We generate \(M_1 = 10\), \(M_2 = 100\), \(M_3 = 1000\), \(M_4 = 10000\) and \(M_5 = 100000\) paths. The paths are plotted in separate figures for the three first cases.

The price dynamics are implemented using the function shown below, which is used to calculate each path. Notice that, since the Wiener process has the distribution \(W_t \sim N(0,t), \forall t \in (0, t]\), we can think of such a process as a cumulative summation of normally distributed random variables. Thus, on our uniform discrete grid

\[ 0 = t_0 < \ldots < t_n = T, \] we can simulate the following relation

\[ X(t_0) = 0, \hspace{0.3em} X(t_1) \sim N(0,t_1) = N(0, t_0+h), \ldots, \hspace{0.3em} X(t_n) \sim N(t_0 + nh) = X(t_1) + \ldots + X(t_{n-1}), \] where \(X\) represents the Brownian motion, \(n\) is the number of steps in the discretization and \(h = \frac{T}{n}\) is the step size.

VIL HAN HA EULER-MARUYAMA HER OG IKKE DENNE FORMELEN? DET ER ENKELT Å BYTTE UT, MEN MÅ SPØRRE (FOR EM ER IMPLEMENTERT LENGER NEDE UANSETT!)

# Price path stochastic process. 
price.path <- function(s0, t, mu, sigma){
  Wt <- rnorm(n = length(t), sd = sqrt(t[2]-t[1])) # Draw length(t) normally distributed variables with mean 0 and standard deviation h. 
  W2 <- cumsum(Wt) # Step size is constant, grid is uniform. We calculate the cumulative sum of the standard normal draw to simulate the Wiener process, which is N(0,t) at time t, i.e. the variance increases when the process is run further and further away from t = 0. Notice that we also multiply by the (uniform) stepsize used on the (uniform) grid, in order to transform the N(0,1) values to N(0,t). 
  s0*exp(mu*t)*exp(sigma*W2-sigma^2/2*t)
}

# VIL HAN EGENTLIG HA EULER-MARUYAMA HER TIL Å LAGE HVER PATH I STEDET!? ELLER ER DETTE OK?

The Monte Carlo estimator for \(\hat{S}_T\) is calculated separately for each of the values of \(M_i, \hspace{0.5em} i \in \{1,2,3,4,5\}\). The 95% confidence interval (CI) is provided for each estimator. A comparison between these estimators and the analytical solution of \(\mathbb{E}(S_T)\) is done and differences are explained.

The Monte Carlo estimator for \(\hat{S}_T\) is simply calculated by averaging the values of all the different BSM price paths plotted above at time \(T = 1\). The \((1-\alpha)\% = (1-0.05)\% = 95\%\) CIs are calculated by finding the standard error of the values of all the different BSM price paths at time \(T\) and using the approximation given by

\[ CI_{\alpha} = \left(\hat{S}_T - z_{\alpha/2}\frac{se_{\hat{S}_T}}{\sqrt{M}}, \hat{S}_T + z_{\alpha/2}\frac{se_{\hat{S}_T}}{\sqrt{M}}\right), \] where \(se_{\hat{S}_T}\) is the aforementioned standard error, \(z_{\alpha/2}\) is the \(1-\alpha\) quantile of the standard normal distribution and \(M\) is the number of price paths simulated. This is an asymptotically valid \(1-\alpha\) CI, which means it becomes closer to the exact analytic value when \(M\) is increased.

Monte Carlo Estimation for S_T, varying M
M1 = 10 M2 = 100 M3 = 1000 M4 = 10000 M5 = 100000
Est. 12.79035 12.27904 12.41963 12.32897 12.36536
Lower CI 11.86917 11.90612 12.28574 12.28727 12.35226
Upper CI 13.71154 12.65195 12.55352 12.37067 12.37845

Now, what is the analytical solution for \(\mathbb{E}(S_T)\)? We know that the process of the risky asset in \(t \in [0,T]\) is distributed as

\[ S_t \sim \mathcal{L}N(\mu^{*}, \sigma^{*2}), \]

where \(\mathbb{E}(S_T) = \mu^{*}\). In addition, we know that, when \(W_t \sim N(0,t)\), it follows that

\[ X_t \sim N(\ln S_0 - \left(\frac{\sigma^2}{2} - \mu\right)t, \sigma^2t) = N(\mu_{X_t}, \sigma^2_{X_t}), \] where \(S_t = e^{X_t}\). Finally, we know that the expected value of the log-normally distributed variable \(S_t\) is \(\exp{\left(\mu_{X_t} + \frac{\sigma^2_{X_t}}{2}\right)}\). This means that the analytical solution for \(\mathbb{E}(S_T)\) is

\[ \mu^* = \exp{\left(\mu_{X_t} + \frac{\sigma^2_{X_t}}{2}\right)} = \exp{\left(\ln S_0 -\left(\frac{\sigma^2}{2} - \mu\right)t + \frac{\sigma^2t}{2}\right)} = S_0 e^{\mu t}, \] which in this case has the numerical value

#> [1] 12.36545

From the results above we can clearly see that the MC estimations move closer and closer to the analytical solution when the number of paths \(M\) is increased. For \(M_5\) the estimation is precise to the first 3 decimals, which I would regard as a very good estimation. We also see that the CI’s clearly change with the increase of \(M\); the lower value of the CI’s increase with \(M\) and the upper value of the CI’s decrease with \(M\), meaning that the \(95\%\) CI’s become narrower when the number of paths increase. This is what we expect from the MC theory, since the estimators are unbiased and, as stated by the strong Law of Large Numbers, the sample average converges a.s. to the true expected value.

Now we fix the number of paths \(M^* = 1000\) and vary the values of \(n\), i.e. the number of steps, while repeating the discussion done above.

Monte Carlo Estimation for S_T, varying n
n1 = 12 n2 = 24 n3 = 250 n4 = 1000
Est. 12.38563 12.27697 12.43257 12.40517
Lower CI 12.24430 12.15002 12.30242 12.27156
Upper CI 12.52696 12.40391 12.56272 12.53879

We can make several similar observations in this case;

Notice however that is is the number of paths \(M\), and not the number of discretization points \(n\), that yields dramatic differences when the value is increased. In other words, the estimates move closer to the true value when \(n\) is increased while \(M\) is fixed, but the changes seem to be more dramatic when \(n\) is fixed and \(M\) is increased. It is however important to notice that, in our experiment, \(n\) is only changed across three orders of magnitude, while \(M\) is changed across five orders of magnitude, which might lead to a somewhat biased observation or discussion.

2 Problem B - Option Pricing I

We calculate the price of the given Call option (parameters are given in code block above) with the Black-Scholes-Merton formula.

#> [1] 3.149899

As we can see, the price is approximately equal to 3.15, when rounded to 2 decimals.

We implement a Monte Carlo estimator for this option price by simulating paths with the Euler-Maruyama method for given steps and paths.

Notice that for the path-independent options, like standard European call and put options, it is not needed to save the price path as is done below. This is done to keep the function as general as possible, in order to re-use it for the rest of the assignment. THIS COULD EASILY BE CHANGED IN THIS CASE IF TOO SLOW! (simply overwriting the previous value in the loop in every iteration).

Assuming that \(\mu = r\).

Relative Errors
n = 10 n = 100 n = 1000
M = 10 -0.5821167 0.5566707 0.6125545
M = 100 0.0255317 0.0021910 -0.3364722
M = 1000 0.0522963 0.0273627 0.0263294
M = 10000 -0.0066211 0.0508549 0.0291815
M = 100000 -0.0035652 -0.0019417 0.0029947
Absolute Errors
n = 10 n = 100 n = 1000
M = 10 1.8336091 1.7534568 1.9294852
M = 100 0.0804222 0.0069013 1.0598536
M = 1000 0.1647280 0.0861897 0.0829349
M = 10000 0.0208557 0.1601878 0.0919187
M = 100000 0.0112299 0.0061162 0.0094331

Interpretation of results in view of \(n\) and \(M\).

TESTER MED EKSEMPLER I SLIDES!

3 Problem C - Option Pricing II

A Monte Carlo estimator is implemented for a specific Asian Call Option. This Asian Call Option averages prices every Friday. We set \(n = 252\) and assume that \(n = 1\) is Monday. This means that we average the prices \(t_5, t_{10}, t_{15}, \ldots\). The paypff profile for this option is thus \((\bar{S}-K)^+\), where \(\bar{S}\) is the arithmetic average over all the prices \(t_i, \hspace{0.2em} i \in \{5, 10, 15, \ldots\}\). We set \(M = 10000\) and use the parameter set \((s_0, T, r, \sigma, K) = (20, 1, 0.02, 0.24, 20)\).

As we can see from the result above, the price estimated via the MC estimator is 1.17. THIS PRICE CAN PERHAPS BE CHECKED WITH A CALCULATOR ONLINE!

A Monte Carlo estimator is implemented for a Lookback option with payoff profile \((S_{max}-K)^+\), where \(S_{max}\) refers to the maximum price during the time to maturity of the option. We use the parameter set \((s_0, T, r, \sigma, K) = (22, 2, 0.02, 0.29, 20)\) THIS K IS NOT SPECIFIED IN THE PROBLEM DESCRIPTION, BUT I WILL JUST LEAVE IT HERE NOW!

As we can see from the result above, the price estimated via the MC estimator is 0. THIS PRICE CAN PERHAPS BE CHECKED WITH A CALCULATOR ONLINE!